Ultracold dipolar gas in an optical lattice: the fate of metastable states 
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We study the physics of ultracold dipolar bosons in optical lattices. We show that dipole-dipole 
interactions lead to the appearance of many insulating metastable states. We study the stability and 
lifetime of these states using a generalization of the instanton theory. We investigate also possibilities 
to prepare, control and manipulate these states using time dependent superlattice modifications and 
modulations. We show that the transfer from one metastable configuration to another necessarily 
occurs via superfluid states, but can be controlled fully on the quantum level. We show how the 
metastable states can be created in the presence of the harmonic trap. Our findings open the way 
toward applications of the metastable states as quantum memories. 
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I. INTRODUCTION 

Ultracold dipolar gases have recently attracted a lot 
of attention, both from the experimental and theoretical 
side [H, [3, H, 0| ■ Experiments show that it is possible to 
cool Chromium below the degeneracy temperature p|, 
and by using a Feshbach resonance one can reduce the 
s-wave contact interaction [6( such that the physics of 
the system is dominated by the anisotropic dipole-dipole 
interaction between atoms. These achievements together 
with the rather remarkable progress in cooling and trap- 
ping of dipolar molecules p} clearly indicates that these 
systems indeed are within experimental reach. Particu- 
larly interesting in this context are ultracold dipolar gases 
in optical lattices [2|, |8( , which offer novel possibilities of 
studying strongly correlated states. 

Theoretical studies 0, [lofl . and experiments [HI on 
bosons in optical lattices have pointed out the existence 
of two main kinds of phases: (i) a superfluid phase (SF), 
characterized by a uniform non-zero order parameter, in 
which particles are delocalized over the whole lattice, and 
(ii) a Mott-insulator phase (MI) of localized atoms, in 
which the order parameter is zero. As soon as one in- 
troduces a long range interaction between atoms, new 
phases appear, both in the SF [H, El and MI [Tj, Gj| 
region of the phase diagram: the supersolid phase (SS), 
which features a non-zero order parameter following a 
modulated pattern, and the charge density wave (CDW) 
in which the localized atoms follow modulated patterns 
with a resulting zero order parameter. These effects are 
particularly strong in dipolar gase s in optical lattices [l5[ . 
Moreover, it has been shown [lj| that the lower tunnel- 
ing region of the phase diagram is characterized by the 
existence of many almost degenerate metastable states of 
Mott-like distribution of atoms in the lattice. These dis- 
tributions consist of localized atoms following a specific 
pattern with a filling factor (average number of atoms 
per site) which is in general not integer. 

The large number of metastable states suggests the 
analogy between dipolar gases and classical complex sys- 
tems such as neural network models [l7| . or spin glass 
models [l8j . As it is well known, classical complex sys- 
tems, and in particular neural networks, may serve very 



well as very efficient classical distributed memory mod- 
els. They are robust with respect to the damage of part 
of the network, and they "recognize patterns with dis- 
tortion" , acting as associative memory. It is natural to 
ask if these properties could not turn out to be useful 
also for quantum memories. Quantum memories serve to 
store in a robust way quantum states, i.e. not only some 
"classical" patterns, but also quantum fluctuations (see 
for instance [19|). It would be interesting to combine the 
best of the two worlds: robustness, associativity and large 
storage capacity of " classical" distributed memories, with 
quantum stability for the storage of fluctuations. This is 
the motivation and far reaching goal of this paper. The 
paper itself concentrates on the first steps toward this 
goal: stability, control, preparation and manipulation of 
the metastable states in ultracold dipolar gases in optical 
lattices. 

We focus on the insulating metastable states of the 
system. We show how the appearance of these states 
crucially depends on the dipole-dipole interaction. Our 
calculations on the stability of these states show that 
their lifetime strongly increases when hopping is sup- 
pressed and scales exponentially with the number of sites 
involved in the tunneling process to other metastable 
states. We show that once the system is prepared in a 
certain metastable configuration, it is necessary to pass 
through the SF region of the phase diagram in order to 
dynamically pass from the given metastable configura- 
tion to another one, and that this is a quantum controlled 
process. 

The paper is organized as follows. In Sect. |TT]we in- 
troduce the model. In Sect. Hill we first derive the mean 
field (MF) Hamiltonian, then calculate the ground and 
metastable states of the system, for the case in which the 
dipole-dipole interaction is a small perturbation with re- 
spect to the contact interaction (U/Unn = 20, see text 
for details). This section contains several subsections: 
we study the behavior of the system with respect to the 
cut-off range of the dipole-dipole interaction and the size 
of the elementary cell that reproduces the infinite lattice. 
Non-uniform lattices are also discussed. Low energy ex- 
citations are discussed in Sect. IIV1 while the stability of 
the metastable states is discussed in Sect. [Vj through an 
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instanton approach. In Sect. IVII we study how to ma- 
nipulate, in a deterministic way, the metastable config- 
urations by changing in time the lattice parameters. In 
Sect. IVII| we treat the effects of a confining trap on the 
system. We discuss our results in Sect. IVIIII 



range interaction to be a small perturbation with respect 
to the contact interaction U/Unn = 20, because this 
regime is reachable with Chromium atoms without the 
need for extreme modification of the scattering length. 



II. THE MODEL 



III. GROUND STATE AND METASTABLE 
STATES 



We study a single component gas of bosons (i.e. spin or 
pseudo-spin, polarized) [ijl, H(| in an optical lattice. We 
assume the temperature of the system to be low enough 
such that we can restrict to the first Bloch band, and the 
system is well described by the extended Bose-Hubbard 
Hamiltonian: 



H = 



-gE ( a * a i + aia l) ~ zZ ^ 

E^ n i( n t - 1) + E E -jniThj, 

« I <(y»«- 



(1) 



where J is the tunneling coefficient, U the on-site interac- 
tion, U jr the strengths of the dipole-dipole interaction at 
different relative distances, and fi the chemical potential 
which fixes the average atomic density. In our notation 
(ij) represents nearest neighbors, and ((ij))jr represents 

neighbors at distance i. 

We describe our system with Hamiltonian (TT|), and a 
Gutzwiller ansatz for the wave function 10] 



(2) 



where \i,n) denotes the Fock state of n atoms at site 
i. In particular the time-dependence of the Gutzwiller 
coefficients /A allows to study the evolution of the state 
in real and imaginary (r = —it) time [TBI . l2l| 



.dfk 



dt 



viV^f^i+tiVW+TfZU + (3) 



(i) 



i(nj - 1) + Ugn^rii - firii 
1 



I n j 



with <pi = ($|ai|3>), <pi = J2(j) z <Pj, n* = (®\4 a i\®)> and 



We consider an infinite two dimensional (2D) square 
lattice. We assume the dipoles to be polarized such that 
the atoms in the lattice experience a repulsive dipole- 
dipole interaction in all directions of the plane. The be- 
havior of the system is determined by the parameters of 
Hamiltonian |T]), in particular the ratio U/Unn between 
the contact interaction and the strength of the first near- 
est neighbor (INN) dipole-dipole interaction. In exper- 
iments with Chromium atoms, by using Feshbach reso- 
nances, it is possible to control the ratio U /Unn and even 
turn it down to zero Q . In this paper we assume the long 



Consider an infinite 2D square lattice reproduced by a 
4x4 elementary cell with periodic boundary conditions, 
filled with dipolar atoms. Such a system is described by 
the Hamiltonian ([1]) and we know that it is characterized 
by the existence of many almost degenerate metastable 
states [161 ]. 

Using imaginary time evolution in Eq. (J3J) , it is pos- 
sible to find the ground state of the system. However 
many times this process gets stucked in local minima of 
energy and in general it is very difficult to reach the ac- 
tual ground state [l6j]. This is a clear signature of the 
existence of metastable states. To find all the metastable 
states, we use a combined mean field and perturbative 
approach. 

We want to write Hamiltonian (TT|) as a sum of single- 
site Hamiltonians. Writing the annihilation operator as 
di = hi + (fi, we can perform the mean field decoupling 
on the product 



a-a,j — a- tp j + ajLpi + ipupj + a\a,j 



(4) 



where in the last step we have assumed small fluctua- 
tions, characteristic of the Mott or the deep superfluid 
states, and replaced a\a,j ~ 0. In Hamiltonian ([T]) we 
now replace ajaj with the expression calculated above, 
and find the mean field Hamiltonian 



Hmf — H Q + Hi , 



where 



fin, + — ni(n t - I) + -^ n i,i n i 



Ho = E 

i 



(5) 

(6) 
(7) 



and we have neglected terms of the order of ipf . 

Given a classical distribution of atoms in the lattice 
1$/), e.g. (I) in Fig. [2 that fulfills ff |$j) = -Eo|*/>, 
we want to know whether this configuration is stable 
or not with respect to particle- hole excitations [§]; this 
holds not only for being the ground state, but also a 
metastable state. For each stable configuration there is a 
region in the fx — J plane, called Mott lobe, in which the 
order parameter is zero due to the perfect localization 
of the atoms at the lattice sites. Therefore, to calcu- 
late the Mott lobe of 1$/), we have to evaluate the order 
parameter tfi = (a*) = Tr(ai/?) at each site of the lat- 
tice. The partition function Z = Ti{er" MF ), after a 
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Dyson expansion of the exponential, is in the lowest rel- 
evant order Z ~ Tr(e — ^°), with (3 being the inverse of 
temperature. The MF density matrix is p = ^e~@ HuF , 
which in the limit of zero temperature (/3 — > oo) can be 
expanded around Eq, and becomes p = e P E o e -0H M F ^ 

Using again the Dyson expansion of the exponential, 
we obtain the order parameter as 



0E O / Tr 



dr 



J(p^ E ° I Tr 
/o 



dr. (8) 



Performing the integral (jHJl in the zero temperature limit, 
we trace around 1$/). To avoid the divergence of the 
integral one has to require that by adding (removing) 
one particle to (from) at any site i, as shown in the 
Fig. [TJ the energy increases, i.e. that the state |<frj) is 
a minimum with respect to particle-hole excitations in 
some range of the parameters J and p. 



(i) 



(la) 



(lb) 



m 



are never satisfied. This is because the requirement that 
by adding (removing) one particle to (from) at any 
site i the energy increases is false, and the integral ([5]) 
indeed diverges. 

We follow the same procedure for every possible classi- 
cal distribution of atoms in the lattice, for filling factors 
v = N a /N s (number of atoms per number of sites) rang- 
ing from v = 1/16 to v = 1. Figs. (a, to, c) show the 
phase diagram calculated in this way for a range of the 
dipole-dipole interaction cut at the first (INN), second 
(2NN) and fourth (4NN) nearest neighbor (see also [III]). 
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FIG. 1: Gray sites are occupied by one atom and white sites 
are empty. A "classical" distribution of atoms in the 

lattice (I). The same distribution with one additional atom 
(la) and one removed atom (lb). 



After simple algebra, one finds the order parameter to 
fulfill 



1 



Un, 



dip 



U(n % - 1) - p + V 



dip 



,(9) 



where V^* is the dipole-dipole interaction of one atom 
placed at site i with the rest of the lattice, and the con- 
ditions for convergence are 



(ni-l)U + V^<n<n i + V^ 



dip ' 



(10) 



One finds such an equation and conditions (flT))) for 
every site i of the lattice. The convergence conditions 
are simple and among them one has to choose the most 
stringent to find the boundary of the lobe at J — 0. In- 
stead the equations for the order parameters are coupled 
due to the (pi term, which can be written in a matrix 
form M(p, U, J) ■ tp — 0, with cp = (■ ■ ■ <p%- ■ ■), and have 
a non trivial solution. For every p, the smallest J for 
which det [M(p, U, J)] = gives the lobe of configura- 
tion |$)/ in the p — J plane. Notice that if the chosen 
configuration is not stable, one finds that conditions (TIT 



FIG. 2: (Color online) (a,b,c) Phase diagram with a range 
of the dipole-dipole interaction cut at the first, second and 
fourth nearest neighbor respectively. The thick line is the 
ground state and the other lobes correspond to the metastable 
states, the same color corresponding to the same filling factor. 
In (c) filling factors range from v = 1/8 to v = 1. Metastable 
configuration appearing at the first nearest neighbor (I), and 
second (Ila-IIb), and the corresponding ground state (GS); 
the metastable states remain stable for all larger ranges of 
the dipole-dipole interaction. 



For 4NN, shown in Fig. [2] (c), the low tunneling region 
of the phase diagram consists of many Mott insulating 
states with different filling factors, ranging from v = 1/8 
to v = 1. The ground state (thick line) is a multiple 
of even filling factors (even number of atoms in the lat- 
tice), while metastable states (thin line) show up also 
with odd values of filling factors apart from 1/16 and 
15/16. The phase diagram presents almost perfect parti- 
cle hole duality induced by the strong value of the on-site 
interaction, meaning that configurations at filling factors 
bigger than v = 1/2 show up with the same number and 
the same distribution of holes in the lattice as the num- 
ber and distribution of atoms for filling factors smaller 
than v = 1/2. Some of the metastable configurations 
at v = 1/2 are presented in Fig. [2] (I,IIa,IIb) with the 
corresponding ground state (GS). We also find supersolid 
domains in the superfluid region of the phase diagram, 
but we do not consider this issue here. 
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FIG. 3: (Color online) Boundaries of the GS Mott lobes at 
zero tunneling, calculated for even filling factors from 2/16 
to 1/2, as a function of the range (RNN) of the dipole-dipole 
interaction. The colors are the same as for the previous figures 
of the lobes. Notice the discontinuity in the GS after RNN = 
4 that will be filled by other fractional filling factors. 

A. Range of interaction 

When the dipole-dipole interaction is absent (Unn = 
0), the phase diagram in the low tunneling region, is given 
by Mott insulator lobes MI(n) with exactly n particles 
per site depending on the value of the chemical potential 
[i @ . In Fig. [2] we have shown the phase diagram for a 
range of dipole-dipole interaction that is cut at the first 
(INN) (a) and second (2NN) (b) nearest neighbor. No- 
tice that as the range of interaction increases, the lower 
point of the MI(l) lobe at J = moves up and lobes for 
fractional filling factors appear in the lower part of the 
phase diagram. The long range interaction is responsible 
for the appearance of the lobes under the MI (I) lobe. 
The checkerboard starts to appear already at INN with 
a small number of metastable insulating lobes and if the 
range of interaction increases the checkerboard moves up 
in the phase diagram as shown in Fig. (Hb). 

As the range of interaction increases new fractional 
filling factors appear. For instance in a 4 x 4 elementary 
cell the smallest allowed filling factor is 1/16. By cutting 
the long range interaction at 4NN we observe that filling 
factor 1/16 is not present, because following conditions 
(fTT)|) . the configuration of one atom in the 4 x 4 elementary 
cell is not stable with respect to particle-hole excitation. 
Nevertheless this configuration becomes stable for larger 
ranges of the dipole-dipole interaction. 



In Fig. [3] we plot the boundaries at zero tunneling of the 
ground state insulating lobes with filling factors multiple 
of 2/16, from 2/16 to 1/2 as a function of the range of the 
interaction. Notice that for RNN < 4 the ground state 
covers entirely the fi domain, from /x = up to the max- 
imum value of filling factor 1/2 zero tunneling boundary, 
while for larger values of the range of interaction discon- 



tinuities start to appear; these are filled by other frac- 
tional filling factors. Notice also that the boundaries at 
zero tunneling stabilize to steady values when the range 
of dipole-dipole interaction is sufficiently large. 



B. Size of the elementary cell 

The size of the elementary cell also plays an important 
role in the allowed filling factors. Indeed, in a N x N cell 
it is impossible to see filling factors smaller than 1/N 2 . 
In order to see the ground state of filling factor 1/2, a 
2x2 elementary cell is sufficient, but to be able to see 
filling factors close to zero and, as a consequence of the 
particle hole duality, close to 1, one has to increase the 
size of the elementary cell. 

Given a range of the dipole-dipole interaction, there is 
a rule of thumb to find which is the smallest GS filling 
factor allowed. It consists of placing atoms in an infinite 
lattice at the smallest possible interatomic distance com- 
patible with zero dipole-dipole interaction in the system, 
and find the dimension of the elementary cell compat- 
ible with this atomic distribution. Table [T] shows the 
relation between the cut-off range of dipole-dipole inter- 
action RNN and the GS minimal filling factor vgs, for 
a cut-off range of interaction up to the eighth nearest 
neighbor. The corresponding lobes in the fi — J plane are 
shown in Fig. [U 
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TABLE I: Range of the dipole-dipole interaction RNN and its 
corresponding GS minimal filling factor vgs- 
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FIG. 4: Lobes of the GS minimal filling factor for RNN = 
1, 8. As the range of dipole-dipole interaction increases the 
tip of the corresponding lobe gets smaller. 
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Notice that as the cut-off range of the long range inter- 
action increases, \x m ax shows a tendency to decrease as 
well as the tip of the lobe for the corresponding minimal 
filling factor. 

In the following, unless differently specified, we will 
consider the elementary cell to be 4 x 4 with periodic 
boundary conditions and the range of dipole-dipole in- 
teraction cut at the fourth nearest neighbor to cover a 
sufficiently round region of interaction in the lattice. 

C. Non-uniform lattices 

Another interesting thing is to see what happens to 
an insulating lobe when we add to the two dimensional 
lattice a superlattice, mimicked by a local chemical po- 
tential Afii with a specific pattern. This can be useful 
for applications such as initialization and manipulation 
of the metastable states. 

We replace the chemical potential in Eqs. (|9l fT0|) with 
fi —¥ fj, — A//j. Our convention is that A/ij < for a 
deeper well, such that it is energetically favorable for an 
atom to stay in it. For any choice of the A/x^, one can 
easily calculate the effect of the superlattice to a given 
insulating lobe. 



towards fi — as Afi decreases. This is easily understood 
by looking at inequalities (jlOp . that in this case become 

(m - i)u + v$ + Afii < m < m + vj£ + A/ii. (ii) 

The upper limit is given by the smallest of the right hand 
side (r.h.s.) in conditions (fTTj) , i.e. at = where 
A/i; = 0, while for the lowest limit, we have to choose 
the biggest of the left hand side (l.h.s.) condition, where 
rii = 1 and A/i,; < 0. 

In Fig. [5] (b), the local chemical potential follows a 
stripe pattern A/i; = A/i n[ S \ where is the density 
distribution of the stripe (S) state (Ha) of Fig. [2j The 
magnitude of A/i is the same as in the above case for the 
thick, dashed and dotted lines. As A/i decreases, the lobe 
becomes smaller due to the distribution of local potential 
energies that do not favor the checkerboard lobe. It is not 
difficult to see that the lower limit, set by the biggest of 
the l.h.s. of conditions ([H}, is given for rii = and 
A/ij = 0, while the upper limit set by the smallest of 
the r.h.s. of conditions (fTl]) . is found where rii = 1 and 
A^ < 0. 



D. 3D lattices 
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FIG. 5: Insulating lobes of the checkerboard state with a 
negative local chemical potential that follow a checkerboard 
pattern (a), and a stripe pattern (b), see text for details. 
The continuous line shows the lobe without any superlattice 
applied. 



In Fig. [5j we show the effect of two different superlat- 
tices on the checkerboard insulating lobe. In Fig. [5] (a) 

A/i; = A/j, n\ , with A/i < and n\ the density of 
the checkerboard. The thick, dashed, and dotted lines are 
for Ah/Unn = 0,-0.6,-1.2 respectively. As A/i grows 
in magnitude, the lobe becomes bigger as expected. No- 
tice that the upper point of the lobes at J — 0, does not 
change while changing A/i, while the lower point moves 



Optical lattices in real experiments are in general three 
dimensional and one should take into account that atoms 
can tunnel in all directions as well as the anisotropic 
dipole-dipole interaction with the whole lattice. While 
it is experimentally feasible to isolate two dimensional 
layers (2D) such that atoms do not tunnel from one layer 
to the neighboring ones, it is not possible to switch off the 
infra-layer dipole-dipole interaction due to its long-range 
character. However, if the direction of the dipoles is per- 
pendicular to the plain of the layers (as in our model), 
the resulting dipole-dipole interaction between different 
plains is attractive. In the Mott phase this makes ener- 
getically favorable to have the same distribution of atoms 
for all layers [13] ■ It would not be difficult to check it by 
using Eq. ([3]) for a three dimensional system and make 
use of the imaginary time evolution technique. This is 
not the purpose of this work and will be done elsewhere. 



IV. LOW ENERGY EXCITATIONS 

The low-lying excitations are creating particles (p) 
and holes (h) in a given metastable configuration. For 
every site i, at J = the excitations are given by 
Ef = Um-n+V£ and = ^-U( ni -1)-V^, where 



n,i is the density at site i. Clearly the hole excitation for 
rii = is unphysical. At finite J, the excitation spectrum 
uj(k) of a metastable configuration, is given by the small 

fluctuations Sfn ^ (t) around the unperturbed metastable 

state coefficients f„ ■ In a Mott state with exactly rrii 
particles at site i, the only non-zero coefficients are given 



by Writing 



fi l) +SfX>(t) in Eq. ©, and 



f(i) 
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taking into account only linear terms in the fluctuations, 
we get 



« 



J 

—ni(rii - 1) + ThV^p - firu - x$ 



(12) 



where ~ E 0)s E„ >7 + T (/^/^r + i^rf*) > 



and x 



|mi(mj — 1) + rriiV^* — \im,i is an extra 



/v „ t ^ / 1 * ' dip i * 

phase that we have introduced to eliminate the rotating 
phase of the fm coefficients. The only non-trivial terms 



in Eq. (|12[) are therefore 



Jy/r 



(13) 



and their complex conjugates. It is convenient to study 
Eq. [13] and their complex conjugates in the Fourier do- 
main with 8$\t) = ^je*'^^!^,*), x« being the 



2D vector pointing at site i. After simple algebra one 
finds the Fourier modes to fulfill 



- j]T . v /m 4 (m J + l)a^ ) ; i (-k,i) + 
ft). L 



ikd 



(14) 



id£> +1 (k,i) - £?a£ +1 (k,t) + 



« 



- J 5] ^(m < + l)(m i + l)o^ 1 (k,t) + 



e ik-d<« f ( 15 ) 



with = {±(d, 0), ±(0, d)} being the vector of nearest 
neighbors in the lattice, and d is the lattice spacing. We 
look for stationary solutions of Eqs. (|14I15[) with the 

ansatz al l) (k,i) = ^(tye - *"^)' + u£ ) (k)e <t "( k )*. For 
every site i of the elementary cell, Eqs. (|14I15[) become 



[S*-«(k)] ^-i(k)-JE W( 
[if + ^k)]^ 1 (k)-J£ - ) 
[i?-a;(k)] u « +1 (k)-JE (i) 

[B* + W (k)] &(-k)-JE 



y/miirrij + l)<^* 1 (-k) + ^rn^mju^^k) 

+ i)K- + i)^+i(-k) + Vh + iH'^-jW 



e ik ' d<3> = 



V(m 4 + i)K- + i) u ^ +1 (k) + VK + iK^f-k) 



e lkd<3> = 
e ik.d<« = o 



(16) 



V /m *( m j + 1 )«m ) fl( k ) + \/™^J U m-l(-k) 



(i)* 



= 0. 



This set of 47V 2 equations can be reduced depending 
on the symmetry of the density distribution, like in the 
case of the checkerboard where only two sites are im- 
portant. Eqs. (|16p can be written in a matrix form, 

AI ( U ] =0, and have non-trivial solution only if 

V v */ 

det [M] = 0. The excitation spectrum is then given 
by the positive solutions of the last equation. We have 
checked that Eqs. (fTBf lead to an excitation spectrum 
that perfectly agrees with the one calculated in [l4[ for 
the checkerboard and the MI{n) states. In Fig. [U(a), we 
show the lowest excitation branch of the four metastable 
configurations of Fig. [2l for fj, = 3.3Unn, J — 0-1Unn 
and k x d = k y d = k/ir, in the first Brillouin zone. The 
thick line is for the (CB) state, the dashed, dash-dotted 
and dotted lines are for (I), (Ha) and (lib) states re- 
spectively. At the boundaries of the insulating lobes the 
excitation spectrum w(k = 0) goes to zero. 

A. Oscillations 

In the real time evolution in Eq. ([3]), at a constant den- 
sity, the chemical potential fj, gives only a phase factor, 



therefore the tunneling coefficient J is the only important 
parameter. Suppose at time t the system is described by 
a Gutzwiller state \<j)t), we define the population of the 
metastable state 10ms) as 

P MS {t) = Vl(0M5|0 t )| 2 , (17) 

the N* h root of the fidelity, where N s is the number of 
sites of the elementary cell. This definition has the ad- 
vantage on the fidelity that it does not depend on the 
number of sites, while the simple fidelity would be one 
if and only if \<pt) = \4>ms), an d otherwise depend on 
the dimension of the cell and tend to zero for an infinite 
number of sites 7V S . 

In a metastable state, atoms are perfectly localized at 
the sites of the lattice, and the system is Mott insulator. 
By adding some "noise" in the GW coefficients we ran- 
domly remove population from occupied sites and move 
it to empty ones, conserving the total number of atoms, 
and the system is superfluid. In general, for a given ini- 
tial condition close to a metastable state (meaning that 
the density follows the distribution of the metastable 
state plus some noise) and the tunneling coefficient J 
smaller than the tip J t i P of the metastable insulating 
lobe, in the real time evolution we observe small oscilla- 
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tions around a local minimum of the energy with a multi- 
component frequency v. The frequency of oscillation v 
depends on the exact initial condition and on the tun- 
neling coefficient. As an example in Fig. [6] (b) we show 
the real time dynamics of the population of metastable 
state (I). The thick and dash-dotted lines are calculated 
for J = O.OAUnn and J = Q.\2Unn respectively, the 
approximate oscillation frequencies, in units of Ti = 1, 
are given by v ~ 27r/150 and v ~ 27r/227 respectively. 
The dashed line is calculated for J = 0.08Unn and at a 
larger value of the initial perturbation, so that its oscil- 
lation frequency is v ~ 27r/297. 
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FIG. 6: (a) Lowest excitation spectrum of metastable state 
(GS) (thick), (I) (dashed), (Ha) (dash-dotted) and (lib) (dot- 
ted) of Fig. [2] calculated for fi = 3.3I7jvjv, and J = 0.1?7jvjv- 
Population (b) of the metastable state (I) during real time 
evolution, the thick and dash-dotted lines are for a small 
perturbation of the metastable state for J = 0.04[/jvjv and 
J = 0.12C/jvjv, while the dashed line corresponds to an initial 
big perturbation and J = O.OSUnn', (c) the Mott insulating 
lobe of the state (I). 



In Fig. [6] (c) we plot the insulating lobe of configuration 
(I), and the round spots are placed in correspondence 
of the values of the parameters for oscillations shown in 
Fig. © (b). 



V. STABILITY OF THE METASTABLE STATES 

In Sect. (|IIip we have studied an infinite 2D lattice re- 
produced by a 4 x 4 elementary cell with periodic bound- 
ary conditions. We have shown that polarized dipolar 
bosons in such a system feature many almost degener- 
ate metastable states that are stable against particle- 
hole excitations. Another clear sign of the existence of 
metastable configurations is the fact that in the imagi- 
nary time evolution is very difficult to reach the ground 
state and often the process gets stucked in local minima 
of energy. Therefore, we can think of these states as local 
minima of a potential where a particle can be trapped for 



a certain time that depends on the barrier that separates 
it from another local minimum. 

Before studying the stability of the metastable states 
we remind the simple case of a particle in a double well 
potential. 



(a) (b) 




FIG. 7: Particle in a double well (a), and the instanton (b). 



Being the particle at time t = in the right well (xo) as 
shown in Fig. [7J (a), the probability at time T for the 
particle to tunnel in the left well can be calculated using 
the propagator in imaginary time and a path integral 
approach [24j . The probability amplitude for the particle 
to tunnel is given by 

(-x \e- TH \x ) = sinh(Tw e~ So ), (18) 

where loq is of the order of the frequency at which the 
particle oscillates around the local minimum xq, and So 
is the minimal action along the stationary path that con- 
nects xq to — xo in the inverted potential of Fig. [7] (b), 
called an instanton. 

The corresponding probability amplitude in real time 
is obtained from (|18[) by analytical continuation just by 
replacing T = iT, and one finds that the particle has 
tunneled completely to the left well after a time given by 

Tcj q = |e s °. (19) 

Therefore the tunneling time is known once we know the 
frequency of small oscillations ujq and the action Sq along 
the stationary path in the inverted potential. 

Regarding the metastable states, the analogy of the 
tunneling of a particle in the double well potential is the 
process in which a metastable state tunnels into its com- 
plementary, in which the role of particles and holes is 
exchanged as shown in Fig. [8] (I) . 



Given a metastable configuration defined by its 
Gutzwiller coefficients {fn^}, it is not straightforward 
to identify the barrier that separates it from its comple- 
mentary. For the Gutzwiller wavefunction we look for 
a simple parametrization that allow us to identify the 
metastable states and parametrize in a simple way the 
process of exchanging atoms with holes and viceversa in 
certain lattice sites. We will describe the process of pass- 
ing from one local minimum to another one using only 
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FIG. 8: (I) Exchanging particles with holes, and (II) process 
where only in a region of the lattice (first and third row from 
top) the exchange of particles with holes takes place. 



one variable and its conjugate momentum. The process 
in which a metastable state tunnels into a state different 
from its complementary, as shown in Fig. [5] (II), car- 
ries the complication that initial and final states are not 
degenerate. 



is the conserved quantity. 

Lagrangian (|23p is an equation in 2N$ independent 
variables and their conjugate momenta , where Ns is the 
number of sites of the lattice. We want to reduce the 
number of independent variables to one. 

Consider a simpler case when we have only two sites 
with one particle in the left well (0=occupied) and we 
want to parametrize the process in which the particle tun- 
nels into the right well (E=empty), with the constraints 

on {x'n , Pn^ ) that the normalization is respected and the 
number of atoms is conserved 



n=0 
E 



£^ )2 - p " )2 ) = L 



i=0 



(25) 



A. Parametrization and Ansatz 



We make the following ansatz 



The imaginary time Lagrangian of a system [25j . de- 
scribed by a quantum state |$), is given by 



($|$\ _ ($1$) 



(20) 



The coefficients of the Gutzwiller wavefunction Q in 
general can be complex numbers. We write them in this 
way 



(21) 



where i„ and p n are real numbers. With the last pre- 
scription, the Lagrangian of the system becomes 



C(x®,pW)=-i P«i«+(*|ff(4°,P^)l*>- (22) 

i,n— 

We write it in its canonical form 
l 



n 1 1 



£( x ^n\ Pri ^) ^ t Pji %{ X n ' ^? 

i,n— 

where 

PjP = -ip£ ) > and 



(23) 



(24) 



(E) 



,(°) 



P 



(O) 



--P, 



p(E) _ p (0) 
JE) _ JO) 



(26) 



(o) _ 



-P 



(O) 



The only independent variable is q. Its conjugate momen- 
tum P is a complicated function given by dC/dq. When 
(q, P) — (0, 0) the atom is in the left well (O) while at 
(q,P) = (V2,0) it is in the right one (E). The instanton 
is then the stationary path that joins those two points in 
phase space. The action So is calculated along this path. 

For more complicated cases, as the processes described 
by Figs. [SI there are more than two sites that ex- 
change particles with holes and viceversa. We consider 
only two independent sites (O) and (E), subject to the 
parametrization (g, P) explained above, where (O) is oc- 
cupied by one atom and in (E) there is a hole. The 
remaining occupied sites (j S {O}), behave as the inde- 
pendent (O), while the empty ones (j S {E}) behave as 
the independent (E). We also have to take into account 
potential sites that do not change, as in the example of 
Fig. [H] (II) . These conditions together with ansatz (|2"5"1) 
and (f2"6")l . enter in in Eq. (f2"3"| as Lagrange multipliers, 
and the Hamiltonian becomes 
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H(q,P)=H(q,P) + \ x 



EE 



i n—0 



,(i) 2 _ p(i) 2 \ _ 2 



(B) _ T (C) 



+ A 2 (p<°> + P) + A 3 (4 £ 
+A 4 (pf > - P<°>) + A 5 - g) + A 6 (if> - P) + 

E [ A °o (4 j) - ?) + - 4 0) ) + (p 0) -p)+ ^ - ^ (0) ~ 



ie{o> 

+ E Ho (4 j) - + Af x - g) + tf (H 

je{E} 



U) 



p 



+ E H o (*b w - 2) + a° ^ + < p (i) + vlA ]) 



ie{o} 



- J] [A},o^ } + Aj,x - 2) + ^o/'.F' + A ( 
ie{i} 



(27) 



where the first two lines are for the independent sites (O) 
and (E), the third and fourth correspond to sites forced 
to behave like (O) or (E), and in the last two lines we 
have taken into account also possible conditions for sites 
that do not change. 



B. Action and barrier 

Given an initial and final distribution of atoms in the 
lattice, as for the examples explained above, the proce- 
dure will be: (i) identify the sites that exchange particles 
with holes and viceversa, sites that do not change, and 
adjust Lagrange multipliers in Eq. (|27|k (ii) calculate the 
stationary path that starts at (q, P) = (0, 0) and, for de- 
generate states, ends at (q, P) = (\/2, 0); (iii) calculate 
the action along the path, given by 



5*n = 



f£(q,P)dT=f £(g,P)$, 

J J path q 



(28) 



with q = dH/dP from Eq. J3ZJ). 

In Figs. [U [TQl we show the stationary path (a) that 
connects the metastable state (I) with its complementary 
(III). At J — > (thick line), the path is bigger than at 
J ^ (dashed line), then the orbit reaches its minimum 
extension in correspondence of the tip of the lobe J. In 
Figs. [H [10] (b), we plot the action per site as a function 
of the tunneling coefficient. The action diverges at J — > 
and it reaches its minimum in correspondence of the tip 
of the metastable lobe. In Figs. [H [10] (c), we calculate 
the barrier as —7i(q, P = 0) for J = 0. 



In Figs. [TJJ [T2] we plot the path, the action and the 
barrier for a metastable state (I) that tunnels into a non- 
degenerate one (III). Notice that in contrast with the 
previous case, the paths start at (q, P) = (0, 0) but end 




FIG. 9: Stationary paths (a), action per site (b), and barrier 
(c). The initial state (I) at (q,P) = (0,0) and final state (III) 
at (q,P) = (a/2,0). The state in (II) is at a middle point 
(<7,P) = (1,0). 



at < q < y/2. This is due to the fact that the final and 
initial states do not have the same energy. 



By comparing the above figures we observe that given 
a metastable state, a longer lifetime corresponds to a 
lower energy barrier. Small energy differences between 
the initial and the final states and large regions of the 
lattice undergoing particle- hole exchange in the tunneling 
process contribute to large energy barriers. Hence, in 
general it is more likely for a given state to tunnel into a 
state deeper in energy, e.g. the ground state, than into its 
complementary, which implies the exchange of particles 
with holes in the whole lattice. 
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FIG. 10: Stationary paths (a), action per site (b), and barrier 
(c). The initial state (I) at (q,P) = (0,0) and final state (III) 
at (q, P) = (V2, 0). The state in (II) is at a middle point 
(«,P) = (1,0). 




FIG. 12: Stationary paths (a), action per site (b), and barrier 
(c). The initial state (I) at (q,P) = (0,0) and final state (III) 
at (q,P) = (\/2, 0). The state in (II) is the decay point at 
(?,P)~ (1.36,0). 




(c) 
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FIG. 11: Stationary paths (a), action per site (b), and barrier 
(c). The initial state (I) at (q,P) = (0,0) and final state (III) 
at (q,P) = (\/2, 0). The state in (II) is the decay point at 
(«,P)~ (0.89,0). 



VI. DYNAMICS 

Once the lattice is prepared in a configuration with 
a certain symmetry, the capability of manipulating the 
configuration is essential in order to use the system as 
a quantum memory. Given an initial metastable state, 
it would be nice to change in time the lattice parame- 
ters such that the system evolves in a deterministic way 
towards another chosen metastable configuration. Since 
there are many metastable states and many parameters, 
we study the problem in a simplified scenario of a finite 
2x2 square lattice, and with a cut-off range of the dipole- 
dipole interaction at the second nearest neighbor (2NN). 



In this section, we study how to dynamically pass from 
a given configuration to another one, with two different 



FH FR FF 



FIG. 13: Disposition of two atoms in the lattice. In MF, CB 
and CB* is the ground state, while S1,..,S4 are degenerate 
metastable states. 



methods: (i) adiabatic passage in which we look for an 
adiabatic transfer of one state into another one, and (ii) 
through MF real-time evolution. 



A. Adiabatic passage 

We study the exact Bose-Hubbard Hamiltonian |T]) for 
the 2x2 lattice mentioned above at filling factor 1/2. 
There are 6 possible ways of placing two atoms in the lat- 
tice and are shown in Fig. I13( which provide the basis for 
the Hilbert subspace of the Bose-Hubbard Hamiltonian. 
Moreover, the MF phase diagram of such a system [28[ 
consists of a 2-times degenerate checkerboard (CB,CB*) 
ground state and a 4-times degenerate metastable state 
in which particles are vertically or horizontally aligned 
in the lattice as a stripe (S1,..,S4) pattern. The CB lobe 
and S lobes are shown in Fig. [3^b). 

The Hamiltonian in this basis is non diagonal because 
of tunneling. By diagonalizing the Hamiltonian, we find 
the ground state 



I' Si 



1-4X 2 



(Itteu) + \1>CB*)) (29) 



to be a symmetric combination of all the states of the 
basis, where a; is a function of the tunneling coefficient. 

Now we want to add to the lattice a superlattice mim- 
icked by a local chemical potential 



fj, — A/i; = fj, — SfJ, nf B — A/z n 



(30) 
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where nf B and nf 1 are respectively the density distribu- 
tions of CB and SI state of Fig. Ql 

To transfer a CB state to the metastable SI, the pro- 
cedure is: (i) prepare the system in CB; one has to break 
the symmetry of the ground state (j2"5)) by applying a local 
chemical potential that privileges CB, in this case <5/i < 
and A/i = in Eq. ([30]) . Notice that S/i from now on will 
be kept constant. Then (ii), apply a second local nega- 
tive chemical potential A/i in the position of atoms in SI 
state. The last process is shown in Fig. Q3J where we 
plot (a) the spectrum of the Bose-Hubbard Hamiltonian 
as a function of A/i. 

Populations 





(b) 


„„l 


n n 




I 



SI S2 CB CB« S3 S4 




FIG. 14: (Color online) (a) Eigenenergies of the Bose- 
Hubbard Hamiltonian, thick and dashed line are the ground 
and first excited state, as a function of A/i. (b) The ground 
state of the system is initially prepared in a CB state (filled 
bar), while the first excited state consists of equally populated 
stripe states (empty bars), then an Sl-type local chemical po- 
tential is applied, (c) The final populations is given by an SI 
ground state (filled bar) and a CB first excited state (empty 
bar), (d) The variation of the tunneling coefficient. 



At A/i = the system is already prepared in the CB state 
for Sfj, — —Unn as shown from the population graph Fig. 
[T3](b). As A/i decreases the population is adiabatically 
transferred in the SI state Fig. [141(c). In correspondence 
to the anticrossing we have increased the tunneling coef- 
ficient J as shown in Fig. [14] (d) , in order to increase the 
magnitude of the gap. For a different adiabatic transfer 
involving other states, one has to go through the above 
steps with the appropriate " distribution" of local chemi- 
cal potential in the lattice. 

The drawback of this approach is that for a finite sys- 
tem there is no SF-MI phase transition [23| , therefore no 
concept of MI lobes. Nevertheless we can identify traces 
of the ground state MI lobe for the filling factor 1/2. We 
calculate the spectrum of Hamiltonian (Q]) for the 2x2 
system introduced above at A/ij = but for all filling 
factors, and in Fig. [T5] we plot the eigenenergies versus 
the chemical potential. The different slopes are for the 
different filling factors v = 1/4, 1/2, and 3/4 respectively. 




FIG. 15: Spectrum of the Bose-Hubbard Hamiltonian for J = 
(a), and J = 0.18Unn (b). At J = the ground state 
coincides with the MF ground state. In (b) the shaded areas 
represent the superfluid region calculated in MF, while the 
round spots are the phase boundaries of filling factor v = 1/2 
calculated with the Bose-Hubbard Hamiltonian. 



At J = Fig. [13 (a) there are 4 degenerate eigenstates 
both for filling factor v = 1/4 and v = 3/4, while for 
v — 1/2 the number of eigenstates is six: a twice degen- 
erate ground state, and an excited state manifold of four 
degenerate states. When J becomes non-zero (b) the de- 
generacy breaks, v = 1/4 and v — 3/4 both split into 
three levels while v — 1/2 splits into four. We identify 
two types of eigenstates: (i) Mott-like states for which 
the eigenenergies do not depend on the tunneling coeffi- 
cient J (continuous lines), and (ii) superfluid- like states 
that change their eigenenergies as J increases (dashed 
lines). At J = 0, as shown in Fig. [15] (a), the ground 
state boundaries of v = 1/2 range from /i = ^/2Umn 
to /i = 4:11 nn and coincide with the boundaries of the 
checkerboard calculated in MF. At J ^ 0, as shown in 
Fig. [15] (b), we estimate the boundaries of the Mott- 
like state (first thick line) as the crossing points with the 
superfluid-like ground states of v = 1/4 and v = 3/4 
(dashed lines), plotted as round spots in the graph. The 
shaded area is the MF superfluid region around filling 
factor 1/2. The boundaries of the checkerboard calcu- 
lated in MF, enclosed in the shaded area, are different 
due to the low accuracy of this method. 



B. MF real time evolution 

For large lattices, it is more reliable to look at the 
dynamics in MF. To pass from one configuration to an- 
other, it turns out to be necessary to go into the super- 
fluid region of the phase diagram. Even if at the MI-SF 
transition it is impossible to be adiabatic because of the 
continuous excitation spectrum of the SF phase [14J , for 
a certain range of lattice parameters the process works. 
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We describe the dynamics through Eq. ([3]) in real time. 
Two things are worth noticing: (i) to have a non trivial 
dynamic, one has to prepare an initial state with a non 
zero superfluid parameter ifi ^ 0, and of course J ^ 0. 
This is because the coupling term in Eq. §5§ is directly 
proportional to the order parameter. And (ii), since the 
number of particles is a constant of the motion that is 
fixed from the initial condition, the chemical potential 
/i gives only a phase factor during the evolution. As a 
consequence of the latter, at a constant integer density 
the only important point in the phase diagram is the tip 
Jup of the insulating lobe indicating the phase transition 
at constant density, such that for J smaller than J t i V the 
system shows small oscillations around a local minimum 
of the energy, as presented in Fig. [6j while for larger 
values of J, one finds deep superfluid oscillations. 

With the population defined as in Eq. (fT7)) . we aim to 
transfer population from a given metastable configura- 
tion to another one with a different symmetry, by chang- 
ing the lattice parameters. In the MF regime, the Mott 
insulator states are exact eigenstates of the MF Hamil- 
tonian. The effect of a non uniform lattice on a given 
configuration, is only to change the size of its insulator 
lobe but, as long as the lobe exists, the configuration re- 
mains stable. This is consistent with the definition of 
metastable states as local minima of energy and is con- 
firmed by the lobes of Fig. \5\ This means that there is no 
coupling between them and consequently, by changing a 
local chemical potential A^ as in Eq. (|30|) . there is no 
possibility of having an anticrossing of the type of Fig. flU 
but only a perfect crossing. The only possibility of trans- 
ferring population from a metastable state to another 
one, is by passing through the superfluid region (SF) of 
the phase diagram and enter a different metastable insu- 
lating lobe. It turns out that this is a quantum controlled 
process which is very much sensitive to the exact initial 
conditions and the way the parameters change in time. 

We specifically study the 2x2 system introduced above 
with periodic boundary conditions and again, we want 
to transfer population from the CB state to SI by ap- 
plying the time dependent local chemical potential A/i^ 
of Eq. (|3"0"]) . in favor of SI. We also want to change the 
tunneling coefficient J in time, so as to exit the CB lobe 
and, through the superfluid region, enter into the SI lobe. 
Ideally we want the population of SI at the end of the 
process to be one, Psi{tfin) = 1, but the actual value 
of Psi{tfin) is very much sensitive on the exact values 
the parameters take during the dynamics. Specifically, 
we change the lattice parameters smoothly in time as 



-Ctanh 



C 



(t - to) 



Ctanh 



-.tr 



C 



(31) 



where C — 1.7Unn and to — 60/Unn (in units of h = 1), 
are kept constant, while a is a free parameter that sets 
the maximum slope for this function, and 

Jm — Jo 



J(A M ) 



2 

tanh [-sm„] 
tanh [s (A^t 



(32) 



- tanh [s (A^i 
rrii)] — tanh [ 



- m )] 
-sm,] - 



2 Jo 
J m -Jo 
2J,„ 



with s — 15/Unn and Jo = 0.02Unn constants, rrii — 
—2.6Unn fixes the superfluid to Mott insulator transition 
point at Afj,i n = —2.57Unn, whereas J m and m are free 
parameters related with the maximum value of tunneling 
coefficient in the superfluid region and the point Aii Q 
where the CB ceases to exist. Together with the inte nsity 
I r of the random noise that fixes the initial condition [29j , 
the space 



{a, Ari a , J m ,I r } ■ 



(33) 



of our control parameters is in total 4-dimensional. We 
require two extra things: (i) an initial non uniform lattice 
that lifts the degeneracy between CB and CB* in favor of 
the checkerboard, whose intensity is fixed at Sfi — —U nn, 
and (ii) we isolate the rows of the lattice from tunneling 
to one another. 




FIG. 16: (a) The pulse of local chemical potential as a func- 
tion of time, (b) The step-like function is the tunneling coef- 
ficient as a function of A/j,, while the thick (dashed) line is the 
tip of SI (CB) insulating lobe, (c) Population inversion, from 
CB to SI at the end of the process. Notice the oscillation of 
populations when passing through the SF region of the phase 
diagram. 



In Fig. [11)] we show the dynamics of the transferring 
process for a — 40 x 10~ 3 U% N , Ait a = -OAWnn, 
J m = 0.66Unn and I r — 4 x 10~ 3 . In Fig. [16] (a) we plot 
the pulse of local chemical potential Afj, as a function of 
time. The smoothed step function in Fig. [11)] (b) shows 
the tunneling coefficient as a function of A/i, while the 
dashed (thick) line is the tip of the CB (SI) insulating 
lobe. These pulses drive the population Fig. [11)] (c) of SI 
(thick line) to a steady value of Psi(tfi n ) = 0.992 at the 
end of the process (tfi n = 120/Unn), while the popula- 
tion of the CB state (dashed line) diminishes consider- 
ably. Notice the oscillation of populations when passing 
through the SF region of the phase diagram. Notice also 
that due to the definition (fl7|) , populations do not have 
to sum up to one. 

Having such a precise control on the parameters (|33|) 
is very challenging from the experimental point of view. 
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Nevertheless, such a process is robust if there is a rea- 
sonable range in which the parameters can vary without 
affecting the final result. The goal is of course the popu- 
lation of SI to be as close as possible to 1 at the end of 
the process. We discretize the space of parameters (|3"3")l 
arbitrarily, and for every value of the parameters simu- 
late the dynamics represented in Fig. 1161 The resulting 
statistics is shown in Fig. [T7] 



space in which one always comes through the 0.98 pop- 
ulation threshold. This is true also for larger elementary 
cells, for which we have checked that the transferring 
process works the same way. This means that experi- 
mentally one has the freedom of setting the control pa- 
rameters such that their small fluctuations do not affect 
the transfer process. This makes the specific process of 
population transferring from CB to SI quite robust. 




FIG. 17: Percentage of realizations (a) terminating with an SI 
population bigger than threshold, versus threshold itself; as 
the threshold increases less realizations satisfy the required 
precision, (b-d) Slices of the discretize space of control pa- 
rameters; the spots are for processes finishing in SI with at 
least 0.98 population. In (b) we fix I r = 10 x 10 -3 and 
Afi = 0.45(7jvjv, in (c) and (d) we fix a = 40 x lO- 3 ^^ 
and A/j, — —0A5Unn respectively, for J m = 0.66I7jvjv. 



VII. TRAP EFFECTS 

So far we have considered an infinite lattice repro- 
duced with a4x4or2x2 elementary cell with periodic 
boundary conditions. In real experiments atoms first are 
trapped in a harmonic trap and then the optical lattice 
is raised. Therefore it is important to understand the 
behavior of these systems in the presence of a confining 
harmonic trap. Here we calculate the ground state of 
a finite 20 x 20 square lattice, where we superimpose a 
harmonic potential mimicked by local chemical potentials 
HH HE H3], without periodic boundary conditions. The 
range of the dipole-dipole interaction is cut at the fourth 
nearest neighbor as before. 




1 10 20 1 10 20 



We have observed that there is a lower limit at J m = 
0.6Unn, below which the transferring process does not 
work. For values of J bigger than this limit, almost all 
the realizations end up in SI state with a final population 
bigger than 0.8, but the exact value depends on the con- 
trol parameters of the single realization. In Fig. [T7] (a), 
we fix J m — 0.66Unn and plot the percentage of the 
dynamics with Psi{tfin) coming through a given thresh- 
old as a function of the population threshold itself. As 
we increase the threshold the number of simulations end- 
ing up in SI with a population that overcomes the given 
threshold decreases, up to no simulations ending up at 
the ideal value Psi{tfin) = 1- This is a clear signature 
of a quantum controlled process. Notice however, that 
about 36% of our simulations terminate with SI being 
populated at 0.98. 

In Fig. [17] (b,c,d) we show slices of the hypercube 
defined by the discretized space of control parameters 
([3"3"]> . where the spots are placed in correspondence of 
the values giving a dynamics with Psiitfin) > 0.98. In 
Fig. QUO) we fix I r = 10 x 10~ 3 and A^ = -QAW NN , 
in Fig. [T7] (c) and (d) we fix a = 40 x IQ~ 3 U% N and 
A/i Q = —0A5Unn respectively, for J m = 0.66Unn- 
There is a closed region in the discretized {a, A[i ,I r } 




1 10 20 1 10 20 



FIG. 18: (Color online) Density p(x,y) and superfluid param- 
eter |yj(a;,i/)| 2 in the harmonic trap. 



The harmonic potential for the system in Fig. [18] is 
V(x, y) = jn— [0 - xo) 2 + (y~ Vo) 2 ] , where (x , y ) is 
the centre of the two dimensional isotropic trap. The pa- 
rameters for the system in Fig. (a,b) are /i/t/jvjv = 2.8, 
J/U NN = 0.26, and K = 107 x lO^s" 1 in units of % = 1. 
There is clearly a region around the center of the trap 
where the density y) follows a checkerboard pattern 
and where the superfluid parameter \ip(x, y)\ 2 is zero (see 
Fig. [18] (a, b)). Notice the supersolid-superfluid area that 
surrounds the Mott insulating phase. In Fig. [T5] (c,d) 
the parameters are h/Unn — 3.3, J/Unn = 0.16, and 
K = 3.1 x 10 _3 s _1 in units of H = 1. The density in the 
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center of the trap (see Fig. [TH](c)) follows the metastable 
state atomic distribution of Fig. [19] (II), with a zero su- 
perfiuid parameter (see Fig. [18] (d)), while in the outer 
region of the trap a SF state is present. In Fig. [19] (a), 
the dashed and thick lines represent the extension of the 
previous trapping potentials, respectively for the (a,b) 
and (c,d) case. 




FIG. 19: (Color online) (a) The thick line is the lobe of the 
ground state (I), while the thin line represents the lobe of 
metastable state (II). Thick and dashed vertical lines are the 
extension of the harmonic potential of Fig. [18] (c,d) and (a,b) 
respectively. 



VIII. CONCLUSION 

We have studied a single component gas of dipolar 
bosons in a two dimensional optical lattice. The atoms 
feature a polarized dipole moment perpendicular to the 
plane of the lattice resulting in a long range interaction 
repulsive in every direction of the plane. The dipole- 
dipole interaction range has been truncated at the fourth 
nearest neighbor, and we have considered 4x4 and 2x2 



unitary cells with periodic boundary conditions. We have 
shown that such a system possesses many almost degen- 
erate metastable states often competing with the ground 
state. 

We have studied the stability of these states and have 
shown that the tunneling time scales exponentially with 
the number of sites of the elementary cell of the corre- 
sponding metastable configurations in the lattice, with a 
factor which depends in a complicated way on the hop- 
ping parameter J, the energy difference between the two 
metastable states, and the number of lattice sites in- 
volved in the tunneling. In a previous work [lril ] , we also 
showed how to identify the state in the lattice through 
noise correlation measurements. 

The mean field theory calculations have shown that, 
once the system is prepared in one of the metastable 
states, it is necessary to go into the superfluid region 
of the phase diagram in order to break the symmetry of 
the prepared state and transfer it to another one. Even 
though this is a quantum controlled process, very much 
sensitive to the exact values of the control parameters 
during the dynamics, we have shown that the process is 
rather robust. 

The capability of initializing, reading and manipulat- 
ing these systems makes dipolar bosons in optical lattice 
very promising for applications in quantum information 
as quantum memories. 
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